# Replication file for quantitative analyses in
# Gaikwad, Lin, and Zucker, "Gender After Genocide: How Violence Shapes Long-Term Political Representation"
#
# *NOTE*: DATA FILES DO NOT CONTAIN NON-PUBLIC VARIABLES. SEE APPENDICES FOR DETAILS.
#
# Software: R version 4.1.2 / RStudio 2022.07.1+554 "Spotted Wakerobin" Release (7872775ebddc40635780ca1ed238934c3345c5de, 2022-07-22) for macOS / macOS 13.0
#
# Updated November 2022
#

# *************** #
# *************** #
###### SETUP ######
# *************** #
# *************** #

pacman::p_load(tidyverse,
               broom,
               AER,
               cobalt,
               data.table,
               haven,
               interflex,
               janitor,
               lfe,
               lmtest,
               magrittr,
               mediation,
               sandwich,
               stargazer)

set.seed(211)

# ************** #
# ************** #
###### DATA ######
# ************** #
# ************** #

data_commune <- read.csv("data_commune_public.csv") # commune-level data
data_commune_by_party <- read.csv("data_commune_by_party_public.csv") # commune-party-level data
data_hh <- read.csv("data_hh_public.csv") # household-level data
data_hh2 <- read.csv("data_hh_public.csv") # household-level data (education)

# ********************* #
# ********************* #
###### DESCRIPTIVE ######
# ********************* #
# ********************* #

# >> Proof of concept: grave measure ####
# >>> Table D-1, Model 1 ####
sex_ratio_on_graves01a <- lm(SEXRATIO65OV ~ graveDistanceCount_mean
                             + as.factor(ZONE_NAME_mode),
                             data = data_commune,
                             na.action = na.omit)
sex_ratio_on_graves01a_coef <- coeftest(sex_ratio_on_graves01a, vcov = vcovHC(sex_ratio_on_graves01a))
sex_ratio_on_graves01a_coef

# >>> Table D-1, Model 2 ####
sex_ratio_on_graves02a <- lm(SEXRATIO65OV ~ graveDistanceCount_mean
                             + mean_elev 
                             + latitude 
                             + dist_to_phnom_penh 
                             + min_dist_to_provincial_capital 
                             + dist_to_road_meters_median 
                             + as.factor(ZONE_NAME_mode),
                             data = data_commune,
                             na.action = na.omit)
sex_ratio_on_graves02a_coef <- coeftest(sex_ratio_on_graves02a, vcov = vcovHC(sex_ratio_on_graves02a))
sex_ratio_on_graves02a_coef

# >>> Table D-1, Model 3 ####
sex_ratio_on_graves01b <- lm(SEXRATIO1564 ~ graveDistanceCount_mean
                             + as.factor(ZONE_NAME_mode),
                             data = data_commune,
                             na.action = na.omit)
sex_ratio_on_graves01b_coef <- coeftest(sex_ratio_on_graves01b, vcov = vcovHC(sex_ratio_on_graves01b))
sex_ratio_on_graves01b_coef

# >>> Table D-1, Model 4 ####
sex_ratio_on_graves02b <- lm(SEXRATIO1564 ~ graveDistanceCount_mean
                             + mean_elev 
                             + latitude 
                             + dist_to_phnom_penh 
                             + min_dist_to_provincial_capital 
                             + dist_to_road_meters_median 
                             + as.factor(ZONE_NAME_mode),
                             data = data_commune,
                             na.action = na.omit)
sex_ratio_on_graves02b_coef <- coeftest(sex_ratio_on_graves02b, vcov = vcovHC(sex_ratio_on_graves02b))
sex_ratio_on_graves02b_coef

# >> Table H-1: Summary statistics ####
stargazer(data_commune[, c("fcouncil_n", "fcouncil_p", "fcand_n", "fcand_p", "FHEADS", "F_LIT15",
                              "mean_elev", "latitude", "dist_to_phnom_penh_km",
                              "min_dist_to_provincial_capital_km", "dist_to_road_meters_median_km",
                              "graveDistanceCount_mean", "graveDistanceCount10km_mean", "graveDistanceCount2.5km_mean")])

# *************** #
# *************** #
###### TESTS ######
# *************** #
# *************** #

# > Political outcomes ~ genocide ####

# ---- Candidates on party lists
# >> Table 1, Model 1 ####
comm_ols_base01a_fcand <- lm(fcand_p ~ graveDistanceCount_mean
                             + as.factor(ZONE_NAME_mode),
                             data = data_commune)
comm_ols_base01a_fcand_coef <- coeftest(comm_ols_base01a_fcand, vcov. = vcovHC(comm_ols_base01a_fcand))

# >> Table 1, Model 2 ####
comm_ols_base01b_fcand <- lm(fcand_p ~ graveDistanceCount_mean
                             + mean_elev 
                             + latitude 
                             + dist_to_phnom_penh 
                             + min_dist_to_provincial_capital 
                             + dist_to_road_meters_median 
                             + as.factor(ZONE_NAME_mode)
                             ,
                             data = data_commune,
                             na.action = na.omit)
comm_ols_base01b_fcand_coef <- coeftest(comm_ols_base01b_fcand, vcov = vcovHC(comm_ols_base01b_fcand))

# ---- Councilors
# >> Table 1, Model 3 ####
comm_ols_base01a_fcouncil <- lm(fcouncil_p ~ graveDistanceCount_mean
                                + as.factor(ZONE_NAME_mode),
                                data = data_commune)
comm_ols_base01a_fcouncil_coef <- coeftest(comm_ols_base01a_fcouncil, vcov. = vcovHC(comm_ols_base01a_fcouncil))

# >> Table 1, Model 4 ####
comm_ols_base01b_fcouncil <- lm(fcouncil_p ~ graveDistanceCount_mean
                                + mean_elev 
                                + latitude 
                                + dist_to_phnom_penh 
                                + min_dist_to_provincial_capital 
                                + dist_to_road_meters_median 
                                + as.factor(ZONE_NAME_mode)
                                ,
                                data = data_commune,
                                na.action = na.omit)
comm_ols_base01b_fcouncil_coef <- coeftest(comm_ols_base01b_fcouncil, vcov = vcovHC(comm_ols_base01b_fcouncil))

# >> 2.5km grave site bandwidth ####

# ---- Candidates on party lists
# >> Table E-1, Left, Model 1 ####
comm_ols_base02a_fcand <- lm(fcand_p ~ graveDistanceCount2.5km_mean
                             + as.factor(ZONE_NAME_mode),
                             data = data_commune)
comm_ols_base02a_fcand_coef <- coeftest(comm_ols_base02a_fcand, vcov. = vcovHC(comm_ols_base02a_fcand))

# >> Table E-1, Left, Model 2 ####
comm_ols_base02b_fcand <- lm(fcand_p ~ graveDistanceCount2.5km_mean
                             + mean_elev 
                             + latitude 
                             + dist_to_phnom_penh 
                             + min_dist_to_provincial_capital 
                             + dist_to_road_meters_median 
                             + as.factor(ZONE_NAME_mode)
                             ,
                             data = data_commune,
                             na.action = na.omit)
comm_ols_base02b_fcand_coef <- coeftest(comm_ols_base02b_fcand, vcov = vcovHC(comm_ols_base02b_fcand))

# ---- Councilors
# >> Table E-1, Left, Model 3 ####
comm_ols_base02a_fcouncil <- lm(fcouncil_p ~ graveDistanceCount2.5km_mean
                                + as.factor(ZONE_NAME_mode),
                                data = data_commune)
comm_ols_base02a_fcouncil_coef <- coeftest(comm_ols_base02a_fcouncil, vcov. = vcovHC(comm_ols_base02a_fcouncil))

# >> Table E-1, Left, Model 4 ####
comm_ols_base02b_fcouncil <- lm(fcouncil_p ~ graveDistanceCount2.5km_mean
                                + mean_elev 
                                + latitude 
                                + dist_to_phnom_penh 
                                + min_dist_to_provincial_capital 
                                + dist_to_road_meters_median 
                                + as.factor(ZONE_NAME_mode)
                                ,
                                data = data_commune,
                                na.action = na.omit)
comm_ols_base02b_fcouncil_coef <- coeftest(comm_ols_base02b_fcouncil, vcov = vcovHC(comm_ols_base02b_fcouncil))

# >> 10km grave site bandwidth ####

# ---- Candidates on party lists
# >> Table E-1, Right, Model 1 ####
comm_ols_base03a_fcand <- lm(fcand_p ~ graveDistanceCount10km_mean
                             + as.factor(ZONE_NAME_mode),
                             data = data_commune)
comm_ols_base03a_fcand_coef <- coeftest(comm_ols_base03a_fcand, vcov. = vcovHC(comm_ols_base03a_fcand))

# >> Table E-1, Right, Model 2 ####
comm_ols_base03b_fcand <- lm(fcand_p ~ graveDistanceCount10km_mean
                             + mean_elev 
                             + latitude 
                             + dist_to_phnom_penh 
                             + min_dist_to_provincial_capital 
                             + dist_to_road_meters_median 
                             + as.factor(ZONE_NAME_mode)
                             ,
                             data = data_commune,
                             na.action = na.omit)
comm_ols_base03b_fcand_coef <- coeftest(comm_ols_base03b_fcand, vcov = vcovHC(comm_ols_base03b_fcand))

# ---- Councilors
# >> Table E-1, Right, Model 3 ####
comm_ols_base03a_fcouncil <- lm(fcouncil_p ~ graveDistanceCount10km_mean
                                + as.factor(ZONE_NAME_mode),
                                data = data_commune)
comm_ols_base03a_fcouncil_coef <- coeftest(comm_ols_base03a_fcouncil, vcov. = vcovHC(comm_ols_base03a_fcouncil))

# >> Table E-1, Right, Model 4 ####
comm_ols_base03b_fcouncil <- lm(fcouncil_p ~ graveDistanceCount10km_mean
                                + mean_elev 
                                + latitude 
                                + dist_to_phnom_penh 
                                + min_dist_to_provincial_capital 
                                + dist_to_road_meters_median 
                                + as.factor(ZONE_NAME_mode)
                                ,
                                data = data_commune,
                                na.action = na.omit)
comm_ols_base03b_fcouncil_coef <- coeftest(comm_ols_base03b_fcouncil, vcov = vcovHC(comm_ols_base03b_fcouncil))

# >> Mass grave size tests ####

# ---- Candidates on party lists
# >> Table J-1, Model 1 ####
comm_ols_base04a_fcand <- lm(fcand_p ~ graveBodiesSum_stnd
                             + as.factor(ZONE_NAME_mode)
                             ,
                             data = data_commune,
                             subset = between(graveBodiesSum, 0, graveBodiesSum975))
comm_ols_base04a_fcand_coef <- coeftest(comm_ols_base04a_fcand, vcov. = vcovHC(comm_ols_base04a_fcand))
comm_ols_base04a_fcand_coef

# >> Table J-1, Model 2 ####
comm_ols_base04b_fcand <- lm(fcand_p ~ graveBodiesSum_stnd
                             + mean_elev 
                             + latitude 
                             + dist_to_phnom_penh 
                             + min_dist_to_provincial_capital 
                             + dist_to_road_meters_median 
                             + as.factor(ZONE_NAME_mode)
                             ,
                             data = data_commune,
                             na.action = na.omit,
                             subset = between(graveBodiesSum, 0, graveBodiesSum975))
comm_ols_base04b_fcand_coef <- coeftest(comm_ols_base04b_fcand, vcov = vcovHC(comm_ols_base04b_fcand))
comm_ols_base04b_fcand_coef

# ---- Councilors
# >> Table J-1, Model 3 ####
comm_ols_base04a_fcouncil <- lm(fcouncil_p ~ graveBodiesSum_stnd
                                + as.factor(ZONE_NAME_mode),
                                data = data_commune,
                                subset = between(graveBodiesSum, 0, graveBodiesSum975))
comm_ols_base04a_fcouncil_coef <- coeftest(comm_ols_base04a_fcouncil, vcov. = vcovHC(comm_ols_base04a_fcouncil))
comm_ols_base04a_fcouncil_coef

# >> Table J-1, Model 4 ####
comm_ols_base04b_fcouncil <- lm(fcouncil_p ~ graveBodiesSum_stnd
                                + mean_elev 
                                + latitude 
                                + dist_to_phnom_penh 
                                + min_dist_to_provincial_capital 
                                + dist_to_road_meters_median 
                                + as.factor(ZONE_NAME_mode)
                                ,
                                data = data_commune,
                                na.action = na.omit,
                                subset = between(graveBodiesSum, 0, graveBodiesSum975))
comm_ols_base04b_fcouncil_coef <- coeftest(comm_ols_base04b_fcouncil, vcov = vcovHC(comm_ols_base04b_fcouncil))
comm_ols_base04b_fcouncil_coef

# >> By party ####

# ------- Adding party-level FEs

# >> Table K-1, Model 1 ####
fcand_by_party01a <- lm(fcand_p2 ~ graveDistanceCount_mean
                        + as.factor(ZONE_NAME_mode)
                        + as.factor(party),
                        data = data_commune_by_party)
fcand_by_party01a_coef <- coeftest(fcand_by_party01a,
                                   vcov. = sandwich::vcovCL(fcand_by_party01a,
                                                            cluster = data_commune_by_party$COMM_CODE))
tidy(fcand_by_party01a_coef)

# >> Table K-1, Model 2 ####
fcand_by_party01b <- lm(fcand_p2 ~ graveDistanceCount_mean
                        + mean_elev 
                        + latitude 
                        + dist_to_phnom_penh 
                        + min_dist_to_provincial_capital 
                        + dist_to_road_meters_median 
                        + as.factor(ZONE_NAME_mode)
                        + as.factor(party)
                        ,
                        data = data_commune_by_party,
                        na.action = na.omit)
fcand_by_party01b_coef <- coeftest(fcand_by_party01b,
                                   vcov. = sandwich::vcovCL(fcand_by_party01b,
                                                            cluster = data_commune_by_party$COMM_CODE))
tidy(fcand_by_party01b_coef)

# >> Comparing communes in neighboring KR zones (W/SW) ####

# ---- Candidates on party lists

# >>> Table 2, Model 1 ####
felm(fcand_p ~ graveDistanceCount_mean:ZONE_NAME_mode | 0 | 0 | 0,
     data = data_commune[data_commune$distance_to_border_km <= 10 &
                              data_commune$ZONE_NAME_mode %in% c("Southwest", "West"), ]) -> border_fcand1
# >>> Table 2, Model 2 ####
felm(fcand_p ~ graveDistanceCount_mean:ZONE_NAME_mode
     + mean_elev 
     + latitude 
     + dist_to_phnom_penh 
     + min_dist_to_provincial_capital 
     + dist_to_road_meters_median
     + I(school_count > 0)
     + I(townhall_count > 0) | 0 | 0 | 0,
     data = data_commune[data_commune$distance_to_border_km <= 10 &
                              data_commune$ZONE_NAME_mode %in% c("Southwest", "West"), ]) -> border_fcand2

# ---- Councilors
# >>> Table 2, Model 3 ####
felm(fcouncil_p ~ graveDistanceCount_mean:ZONE_NAME_mode | 0 | 0 | 0,
     data = data_commune[data_commune$distance_to_border_km <= 10 &
                              data_commune$ZONE_NAME_mode %in% c("Southwest", "West"), ]) -> border_fcouncil1
# >>> Table 2, Model 4 ####
felm(fcouncil_p ~ graveDistanceCount_mean:ZONE_NAME_mode
     + mean_elev 
     + latitude 
     + dist_to_phnom_penh 
     + min_dist_to_provincial_capital 
     + dist_to_road_meters_median
     + I(school_count > 0)
     + I(townhall_count > 0) | 0 | 0 | 0,
     data = data_commune[data_commune$distance_to_border_km <= 10 &
                              data_commune$ZONE_NAME_mode %in% c("Southwest", "West"), ]) -> border_fcouncil2

stargazer(border_fcand1,
          border_fcand2,
          border_fcouncil1,
          border_fcouncil2)

# > FHH, literacy ~ genocide ####

# ---- Density of female-headed households
# >> Table 3, Model 1 ####
comm_ols_base01a_fheads <- lm(FHEADS ~ graveDistanceCount_mean
                              + as.factor(ZONE_NAME_mode),
                              data = data_commune)
comm_ols_base01a_fheads_coef <- coeftest(comm_ols_base01a_fheads, vcov. = vcovHC(comm_ols_base01a_fheads))

# >> Table 3, Model 2 ####
comm_ols_base01b_fheads <- lm(FHEADS ~ graveDistanceCount_mean
                              + mean_elev 
                              + latitude 
                              + dist_to_phnom_penh 
                              + min_dist_to_provincial_capital 
                              + dist_to_road_meters_median 
                              + as.factor(ZONE_NAME_mode)
                              ,
                              data = data_commune,
                              na.action = na.omit)
comm_ols_base01b_fheads_coef <- coeftest(comm_ols_base01b_fheads, vcov = vcovHC(comm_ols_base01b_fheads))

# ---- Female literacy rates
# >> Table 3, Model 3 ####
comm_ols_base01a_flit <- lm(F_LIT15 ~ graveDistanceCount_mean
                            + as.factor(ZONE_NAME_mode),
                            data = data_commune)
comm_ols_base01a_flit_coef <- coeftest(comm_ols_base01a_flit, vcov. = vcovHC(comm_ols_base01a_flit))

# >> Table 3, Model 4 ####
comm_ols_base01b_flit <- lm(F_LIT15 ~ graveDistanceCount_mean
                            + mean_elev 
                            + latitude 
                            + dist_to_phnom_penh 
                            + min_dist_to_provincial_capital 
                            + dist_to_road_meters_median 
                            + as.factor(ZONE_NAME_mode)
                            ,
                            data = data_commune,
                            na.action = na.omit)
comm_ols_base01b_flit_coef <- coeftest(comm_ols_base01b_flit, vcov = vcovHC(comm_ols_base01b_flit))

# >> 2.5km bandwidth ####

# ---- Density of female-headed households
# >>> Table E-2, Left, Model 1 ####
comm_ols_base02a_fheads <- lm(FHEADS ~ graveDistanceCount2.5km_mean
                              + as.factor(ZONE_NAME_mode),
                              data = data_commune)
comm_ols_base02a_fheads_coef <- coeftest(comm_ols_base02a_fheads, vcov. = vcovHC(comm_ols_base02a_fheads))

# >>> Table E-2, Left, Model 2 ####
comm_ols_base02b_fheads <- lm(FHEADS ~ graveDistanceCount2.5km_mean
                              + mean_elev 
                              + latitude 
                              + dist_to_phnom_penh 
                              + min_dist_to_provincial_capital 
                              + dist_to_road_meters_median 
                              + as.factor(ZONE_NAME_mode)
                              ,
                              data = data_commune,
                              na.action = na.omit)
comm_ols_base02b_fheads_coef <- coeftest(comm_ols_base02b_fheads, vcov = vcovHC(comm_ols_base02b_fheads))

# ---- Female literacy rates
# >>> Table E-2, Left, Model 3 ####
comm_ols_base02a_flit <- lm(F_LIT15 ~ graveDistanceCount2.5km_mean
                            + as.factor(ZONE_NAME_mode),
                            data = data_commune)
comm_ols_base02a_flit_coef <- coeftest(comm_ols_base02a_flit, vcov. = vcovHC(comm_ols_base02a_flit))

# >>> Table E-2, Left, Model 4 ####
comm_ols_base02b_flit <- lm(F_LIT15 ~ graveDistanceCount2.5km_mean
                            + mean_elev 
                            + latitude 
                            + dist_to_phnom_penh 
                            + min_dist_to_provincial_capital 
                            + dist_to_road_meters_median 
                            + as.factor(ZONE_NAME_mode)
                            ,
                            data = data_commune,
                            na.action = na.omit)
comm_ols_base02b_flit_coef <- coeftest(comm_ols_base02b_flit, vcov = vcovHC(comm_ols_base02b_flit))

# >> 10km bandwidth ####

# ---- Density of female-headed households
# >>> Table E-2, Right, Model 1 ####
comm_ols_base03a_fheads <- lm(FHEADS ~ graveDistanceCount10km_mean
                              + as.factor(ZONE_NAME_mode),
                              data = data_commune)
comm_ols_base03a_fheads_coef <- coeftest(comm_ols_base03a_fheads, vcov. = vcovHC(comm_ols_base03a_fheads))

# >>> Table E-2, Right, Model 2 ####
comm_ols_base03b_fheads <- lm(FHEADS ~ graveDistanceCount10km_mean
                              + mean_elev 
                              + latitude 
                              + dist_to_phnom_penh 
                              + min_dist_to_provincial_capital 
                              + dist_to_road_meters_median 
                              + as.factor(ZONE_NAME_mode)
                              ,
                              data = data_commune,
                              na.action = na.omit)
comm_ols_base03b_fheads_coef <- coeftest(comm_ols_base03b_fheads, vcov = vcovHC(comm_ols_base03b_fheads))

# ---- Female literacy rates
# >>> Table E-2, Right, Model 3 ####
comm_ols_base03a_flit <- lm(F_LIT15 ~ graveDistanceCount10km_mean
                            + as.factor(ZONE_NAME_mode),
                            data = data_commune)
comm_ols_base03a_flit_coef <- coeftest(comm_ols_base03a_flit, vcov. = vcovHC(comm_ols_base03a_flit))

# >>> Table E-2, Right, Model 4 ####
comm_ols_base03b_flit <- lm(F_LIT15 ~ graveDistanceCount10km_mean
                            + mean_elev 
                            + latitude 
                            + dist_to_phnom_penh 
                            + min_dist_to_provincial_capital 
                            + dist_to_road_meters_median 
                            + as.factor(ZONE_NAME_mode)
                            ,
                            data = data_commune,
                            na.action = na.omit)
comm_ols_base03b_flit_coef <- coeftest(comm_ols_base03b_flit, vcov = vcovHC(comm_ols_base03b_flit))

# >> Comparing communes in neighboring KR zones ####

# ---- FHEADS
# >>> Table L-1, Model 1 ####
felm(FHEADS ~ graveDistanceCount_mean:ZONE_NAME_mode | 0 | 0 | 0,
     data = data_commune[data_commune$distance_to_border_km <= 10 &
                              data_commune$ZONE_NAME_mode %in% c("Southwest", "West"), ]) %>%
  tidy()
# >>> Table L-1, Model 2 ####
felm(FHEADS ~ graveDistanceCount_mean:ZONE_NAME_mode
     + mean_elev 
     + latitude 
     + dist_to_phnom_penh 
     + min_dist_to_provincial_capital 
     + dist_to_road_meters_median
     + I(school_count > 0)
     + I(townhall_count > 0) | 0 | 0 | 0,
     data = data_commune[data_commune$distance_to_border_km <= 10 &
                              data_commune$ZONE_NAME_mode %in% c("Southwest", "West"), ]) %>%
  tidy()

# ---- F_LIT15
# >>> Table L-1, Model 3 ####
felm(F_LIT15 ~ graveDistanceCount_mean:ZONE_NAME_mode | 0 | 0 | 0,
     data = data_commune[data_commune$distance_to_border_km <= 10 &
                              data_commune$ZONE_NAME_mode %in% c("Southwest", "West"), ]) %>%
  tidy()
# >>> Table L-1, Model 4 ####
felm(F_LIT15 ~ graveDistanceCount_mean:ZONE_NAME_mode
     + mean_elev 
     + latitude 
     + dist_to_phnom_penh 
     + min_dist_to_provincial_capital 
     + dist_to_road_meters_median
     + I(school_count > 0)
     + I(townhall_count > 0) | 0 | 0 | 0,
     data = data_commune[data_commune$distance_to_border_km <= 10 &
                              data_commune$ZONE_NAME_mode %in% c("Southwest", "West"), ]) %>%
  tidy()

# > Business management, educational attainment ~ genocide ####

# ---- Business management
# >> Table 4, Model 1 ####
biz.fhh.biv <- felm(fbusinessmanager ~ graveDistanceCount | ZONE_NAME | 0 | villageid, 
                    data = data_hh, 
                    subset = fhh == 1,
                    weights = data_hh$hhweight[data_hh$fhh == 1])
# >> Table 4, Model 2 ####
biz.fhh.cov <- felm(fbusinessmanager ~ graveDistanceCount + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid, 
                    data = data_hh, 
                    subset = fhh == 1,
                    weights = data_hh$hhweight[data_hh$fhh == 1])
# >> Table 4, Model 3 ####
biz.mhh.biv <- felm(fbusinessmanager ~ graveDistanceCount | ZONE_NAME | 0 | villageid, 
                    data = data_hh, 
                    subset = fhh == 0,
                    weights = data_hh$hhweight[data_hh$fhh == 0])
# >> Table 4, Model 4 ####
biz.mhh.cov <- felm(fbusinessmanager ~ graveDistanceCount + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid, 
                    data = data_hh, 
                    subset = fhh == 0,
                    weights = data_hh$hhweight[data_hh$fhh == 0])

# ---- Educational attainment
# >> Table 4, Model 5 ####
edu.fhh.biv <- felm(fschyrs ~ graveDistanceCount | ZONE_NAME | 0 | villageid, 
                    data = data_hh, 
                    subset = fhh == 1,
                    weights = data_hh$hhweight[data_hh$fhh == 1])
# >> Table 4, Model 6 ####
edu.fhh.cov <- felm(fschyrs ~ graveDistanceCount + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid, 
                    data = data_hh, 
                    subset = fhh == 1,
                    weights = data_hh$hhweight[data_hh$fhh == 1])
# >> Table 4, Model 7 ####
edu.mhh.biv <- felm(fschyrs ~ graveDistanceCount | ZONE_NAME | 0 | villageid, 
                    data = data_hh, 
                    subset = fhh == 0,
                    weights = data_hh$hhweight[data_hh$fhh == 0])
# >> Table 4, Model 8 ####
edu.mhh.cov <- felm(fschyrs ~ graveDistanceCount + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid, 
                    data = data_hh, 
                    subset = fhh == 0,
                    weights = data_hh$hhweight[data_hh$fhh == 0])

# >> 2.5km bandwidth ####

# ---- Business management
# >> Table E-3, Top, Model 1 ####
biz.fhh.biv_2.5km <- felm(fbusinessmanager ~ graveDistanceCount2.5km | ZONE_NAME | 0 | villageid, 
                          data = data_hh, 
                          subset = fhh == 1,
                          weights = data_hh$hhweight[data_hh$fhh == 1])
# >> Table E-3, Top, Model 2 ####
biz.fhh.cov_2.5km <- felm(fbusinessmanager ~ graveDistanceCount2.5km + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid, 
                          data = data_hh, 
                          subset = fhh == 1,
                          weights = data_hh$hhweight[data_hh$fhh == 1])
# >> Table E-3, Top, Model 3 ####
biz.mhh.biv_2.5km <- felm(fbusinessmanager ~ graveDistanceCount2.5km | ZONE_NAME | 0 | villageid, 
                          data = data_hh, 
                          subset = fhh == 0,
                          weights = data_hh$hhweight[data_hh$fhh == 0])
# >> Table E-3, Top, Model 4 ####
biz.mhh.cov_2.5km <- felm(fbusinessmanager ~ graveDistanceCount2.5km + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid, 
                          data = data_hh, 
                          subset = fhh == 0,
                          weights = data_hh$hhweight[data_hh$fhh == 0])

# ---- Educational attainment
# >> Table E-3, Top, Model 5 ####
edu.fhh.biv_2.5km <- felm(fschyrs ~ graveDistanceCount2.5km | ZONE_NAME | 0 | villageid, 
                          data = data_hh, 
                          subset = fhh == 1,
                          weights = data_hh$hhweight[data_hh$fhh == 1])
# >> Table E-3, Top, Model 6 ####
edu.fhh.cov_2.5km <- felm(fschyrs ~ graveDistanceCount2.5km + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid, 
                          data = data_hh, 
                          subset = fhh == 1,
                          weights = data_hh$hhweight[data_hh$fhh == 1])
# >> Table E-3, Top, Model 7 ####
edu.mhh.biv_2.5km <- felm(fschyrs ~ graveDistanceCount2.5km | ZONE_NAME | 0 | villageid, 
                          data = data_hh, 
                          subset = fhh == 0,
                          weights = data_hh$hhweight[data_hh$fhh == 0])
# >> Table E-3, Top, Model 8 ####
edu.mhh.cov_2.5km <- felm(fschyrs ~ graveDistanceCount2.5km + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid, 
                          data = data_hh, 
                          subset = fhh == 0,
                          weights = data_hh$hhweight[data_hh$fhh == 0])

# >> 10km bandwidth ####

# ---- Business management
# >> Table E-3, Bottom, Model 1 ####
biz.fhh.biv_10km <- felm(fbusinessmanager ~ graveDistanceCount10km | ZONE_NAME | 0 | villageid, 
                         data = data_hh, 
                         subset = fhh == 1,
                         weights = data_hh$hhweight[data_hh$fhh == 1])
# >> Table E-3, Bottom, Model 2 ####
biz.fhh.cov_10km <- felm(fbusinessmanager ~ graveDistanceCount10km + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid, 
                         data = data_hh, 
                         subset = fhh == 1,
                         weights = data_hh$hhweight[data_hh$fhh == 1])
# >> Table E-3, Bottom, Model 3 ####
biz.mhh.biv_10km <- felm(fbusinessmanager ~ graveDistanceCount10km | ZONE_NAME | 0 | villageid, 
                         data = data_hh, 
                         subset = fhh == 0,
                         weights = data_hh$hhweight[data_hh$fhh == 0])
# >> Table E-3, Bottom, Model 4 ####
biz.mhh.cov_10km <- felm(fbusinessmanager ~ graveDistanceCount10km + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid, 
                         data = data_hh, 
                         subset = fhh == 0,
                         weights = data_hh$hhweight[data_hh$fhh == 0])

# ---- Educational attainment
# >> Table E-3, Bottom, Model 5 ####
edu.fhh.biv_10km <- felm(fschyrs ~ graveDistanceCount10km | ZONE_NAME | 0 | villageid, 
                         data = data_hh, 
                         subset = fhh == 1,
                         weights = data_hh$hhweight[data_hh$fhh == 1])
# >> Table E-3, Bottom, Model 6 ####
edu.fhh.cov_10km <- felm(fschyrs ~ graveDistanceCount10km + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid, 
                         data = data_hh, 
                         subset = fhh == 1,
                         weights = data_hh$hhweight[data_hh$fhh == 1])
# >> Table E-3, Bottom, Model 7 ####
edu.mhh.biv_10km <- felm(fschyrs ~ graveDistanceCount10km | ZONE_NAME | 0 | villageid, 
                         data = data_hh, 
                         subset = fhh == 0,
                         weights = data_hh$hhweight[data_hh$fhh == 0])
# >> Table E-3, Bottom, Model 8 ####
edu.mhh.cov_10km <- felm(fschyrs ~ graveDistanceCount10km + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid, 
                         data = data_hh, 
                         subset = fhh == 0,
                         weights = data_hh$hhweight[data_hh$fhh == 0])

# > Education spending ~ genocide ####

# ---- % spending on education
# >> Table 5, Model 1 ####
felm(I(education.spending/total.spending) ~ graveDistanceCount + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid,
     data = data_hh2,
     weights = data_hh2$hhweight) -> reg_cses9_education.spending
# >> Table 5, Model 2 ####
felm(I(education.spending/total.spending) ~ graveDistanceCount + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid,
     data = data_hh2,
     subset = FHH == 1,
     weights = data_hh2$hhweight[data_hh2$FHH == 1]) -> reg_cses9_education.spending_FHH
# >> Table 5, Model 3 ####
felm(I(education.spending/total.spending) ~ graveDistanceCount + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid,
     data = data_hh2,
     subset = FHH == 0,
     weights = data_hh2$hhweight[data_hh2$FHH == 0]) -> reg_cses9_education.spending_MHH

# ---- % spending on female education
# >> Table 5, Model 4 ####
felm(I(education_spending_female/total.spending) ~ graveDistanceCount + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid,
     data = data_hh2,
     weights = data_hh2$hhweight) -> reg_cses9_education_spending_female2
# >> Table 5, Model 5 ####
felm(I(education_spending_female/total.spending) ~ graveDistanceCount + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid,
     data = data_hh2,
     subset = FHH == 1,
     weights = data_hh2$hhweight[data_hh2$FHH == 1]) -> reg_cses9_education_spending_female2_FHH
# >> Table 5, Model 6 ####
felm(I(education_spending_female/total.spending) ~ graveDistanceCount + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid,
     data = data_hh2,
     subset = FHH == 0,
     weights = data_hh2$hhweight[data_hh2$FHH == 0]) -> reg_cses9_education_spending_female2_MHH

# ---- Education spending per woman/girl
# >> Table 5, Model 7 ####
felm(I((education_spending_female/females_positive_spending)/1000) ~ graveDistanceCount + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid,
     data = data_hh2,
     weights = data_hh2$hhweight) -> reg_cses9_education_spending_female
# >> Table 5, Model 8 ####
felm(I((education_spending_female/females_positive_spending)/1000) ~ graveDistanceCount + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid,
     data = data_hh2,
     subset = FHH == 1,
     weights = data_hh2$hhweight[data_hh2$FHH == 1]) -> reg_cses9_education_spending_female_FHH
# >> Table 5, Model 9 ####
felm(I((education_spending_female/females_positive_spending)/1000) ~ graveDistanceCount + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid,
     data = data_hh2,
     subset = FHH == 0,
     weights = data_hh2$hhweight[data_hh2$FHH == 0]) -> reg_cses9_education_spending_female_MHH

# >> 2.5km bandwidth ####

# ---- % spending on education
# >>> Table E-4, Top, Model 1 ####
felm(I(education.spending/total.spending) ~ graveDistanceCount2.5km + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid,
     data = data_hh2,
     weights = data_hh2$hhweight) -> reg_cses9_education.spending_2.5km
# >>> Table E-4, Top, Model 2 ####
felm(I(education.spending/total.spending) ~ graveDistanceCount2.5km + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid,
     data = data_hh2,
     subset = FHH == 1,
     weights = data_hh2$hhweight[data_hh2$FHH == 1]) -> reg_cses9_education.spending_FHH_2.5km
# >>> Table E-4, Top, Model 3 ####
felm(I(education.spending/total.spending) ~ graveDistanceCount2.5km + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid,
     data = data_hh2,
     subset = FHH == 0,
     weights = data_hh2$hhweight[data_hh2$FHH == 0]) -> reg_cses9_education.spending_MHH_2.5km

# ---- % spending on female education
# >>> Table E-4, Top, Model 4 ####
felm(I(education_spending_female/total.spending) ~ graveDistanceCount2.5km + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid,
     data = data_hh2,
     weights = data_hh2$hhweight) -> reg_cses9_education_spending_female2_2.5km
# >>> Table E-4, Top, Model 5 ####
felm(I(education_spending_female/total.spending) ~ graveDistanceCount2.5km + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid,
     data = data_hh2,
     subset = FHH == 1,
     weights = data_hh2$hhweight[data_hh2$FHH == 1]) -> reg_cses9_education_spending_female2_FHH_2.5km
# >>> Table E-4, Top, Model 6 ####
felm(I(education_spending_female/total.spending) ~ graveDistanceCount2.5km + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid,
     data = data_hh2,
     subset = FHH == 0,
     weights = data_hh2$hhweight[data_hh2$FHH == 0]) -> reg_cses9_education_spending_female2_MHH_2.5km

# ---- Education spending per woman/girl
# >>> Table E-4, Top, Model 7 ####
felm(I((education_spending_female/females_positive_spending)/1000) ~ graveDistanceCount2.5km + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid,
     data = data_hh2,
     weights = data_hh2$hhweight) -> reg_cses9_education_spending_female_2.5km
# >>> Table E-4, Top, Model 8 ####
felm(I((education_spending_female/females_positive_spending)/1000) ~ graveDistanceCount2.5km + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid,
     data = data_hh2,
     subset = FHH == 1,
     weights = data_hh2$hhweight[data_hh2$FHH == 1]) -> reg_cses9_education_spending_female_FHH_2.5km
# >>> Table E-4, Top, Model 9 ####
felm(I((education_spending_female/females_positive_spending)/1000) ~ graveDistanceCount2.5km + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid,
     data = data_hh2,
     subset = FHH == 0,
     weights = data_hh2$hhweight[data_hh2$FHH == 0]) -> reg_cses9_education_spending_female_MHH_2.5km

# >> 10km bandwidth ####

# ---- % spending on education
# >>> Table E-4, Bottom, Model 1 ####
felm(I(education.spending/total.spending) ~ graveDistanceCount10km + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid,
     data = data_hh2,
     weights = data_hh2$hhweight) -> reg_cses9_education.spending_10km
# >>> Table E-4, Bottom, Model 2 ####
felm(I(education.spending/total.spending) ~ graveDistanceCount10km + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid,
     data = data_hh2,
     subset = FHH == 1,
     weights = data_hh2$hhweight[data_hh2$FHH == 1]) -> reg_cses9_education.spending_FHH_10km
# >>> Table E-4, Bottom, Model 3 ####
felm(I(education.spending/total.spending) ~ graveDistanceCount10km + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid,
     data = data_hh2,
     subset = FHH == 0,
     weights = data_hh2$hhweight[data_hh2$FHH == 0]) -> reg_cses9_education.spending_MHH_10km

# ---- % spending on female education
# >>> Table E-4, Bottom, Model 4 ####
felm(I(education_spending_female/total.spending) ~ graveDistanceCount10km + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid,
     data = data_hh2,
     weights = data_hh2$hhweight) -> reg_cses9_education_spending_female2_10km
# >>> Table E-4, Bottom, Model 5 ####
felm(I(education_spending_female/total.spending) ~ graveDistanceCount10km + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid,
     data = data_hh2,
     subset = FHH == 1,
     weights = data_hh2$hhweight[data_hh2$FHH == 1]) -> reg_cses9_education_spending_female2_FHH_10km
# >>> Table E-4, Bottom, Model 6 ####
felm(I(education_spending_female/total.spending) ~ graveDistanceCount10km + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid,
     data = data_hh2,
     subset = FHH == 0,
     weights = data_hh2$hhweight[data_hh2$FHH == 0]) -> reg_cses9_education_spending_female2_MHH_10km

# ---- Education spending per woman/girl
# >>> Table E-4, Bottom, Model 7 ####
felm(I((education_spending_female/females_positive_spending)/1000) ~ graveDistanceCount10km + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid,
     data = data_hh2,
     weights = data_hh2$hhweight) -> reg_cses9_education_spending_female_10km
# >>> Table E-4, Bottom, Model 8 ####
felm(I((education_spending_female/females_positive_spending)/1000) ~ graveDistanceCount10km + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid,
     data = data_hh2,
     subset = FHH == 1,
     weights = data_hh2$hhweight[data_hh2$FHH == 1]) -> reg_cses9_education_spending_female_FHH_10km
# >>> Table E-4, Bottom, Model 9 ####
felm(I((education_spending_female/females_positive_spending)/1000) ~ graveDistanceCount10km + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | ZONE_NAME | 0 | villageid,
     data = data_hh2,
     subset = FHH == 0,
     weights = data_hh2$hhweight[data_hh2$FHH == 0]) -> reg_cses9_education_spending_female_MHH_10km

# > Political outcomes ~ genocide*socioeconomic empowerment ####
# >> Figure M-1 ####
interflex(data = as.data.table(data_commune),
          Y = "fcand_p",
          D = "graveDistanceCount_mean",
          X = "FHEADS",
          Z = c("mean_elev", "latitude", "dist_to_phnom_penh_km", 
                "min_dist_to_provincial_capital_km", "dist_to_road_meters_median_km"),
          FE = "ZONE_NAME_mode",
          full.moderate = F,
          estimator = "binning",
          nbins = 3,
          na.rm = T,
          theme.bw = T,
          xlab = "Density of female-headed households",
          ylab = "Marginal effect of grave proximity on women as % party list",
          bin.labs = FALSE) -> bin_cand_covar

interflex(data = as.data.table(data_commune),
          Y = "fcouncil_p",
          D = "graveDistanceCount_mean",
          X = "FHEADS",
          Z = c("mean_elev", "latitude", "dist_to_phnom_penh_km", 
                "min_dist_to_provincial_capital_km", "dist_to_road_meters_median_km"),
          FE = "ZONE_NAME_mode",
          full.moderate = F,
          estimator = "binning",
          nbins = 3,
          na.rm = T,
          theme.bw = T,
          xlab = "Density of female-headed households",
          ylab = "Marginal effect of grave proximity on women as % council",
          bin.labs = FALSE) -> bin_council_covar

# >> Figure M-2 ####
interflex(data = as.data.table(data_commune),
          Y = "fcand_p",
          D = "graveDistanceCount_mean",
          X = "F_LIT15",
          Z = c("mean_elev", "latitude", "dist_to_phnom_penh_km", 
                "min_dist_to_provincial_capital_km", "dist_to_road_meters_median_km"),
          FE = "ZONE_NAME_mode",
          full.moderate = F,
          estimator = "binning",
          nbins = 3,
          na.rm = T,
          theme.bw = T,
          xlab = "Female literacy rate",
          ylab = "Marginal effect of grave proximity on women as % party list",
          bin.labs = FALSE) -> bin_cand_covar_LIT # interaction with covariates

interflex(data = as.data.table(data_commune),
          Y = "fcouncil_p",
          D = "graveDistanceCount_mean",
          X = "F_LIT15",
          Z = c("mean_elev", "latitude", "dist_to_phnom_penh_km", 
                "min_dist_to_provincial_capital_km", "dist_to_road_meters_median_km"),
          FE = "ZONE_NAME_mode",
          full.moderate = F,
          estimator = "binning",
          nbins = 3,
          na.rm = T,
          theme.bw = T,
          xlab = "Female literacy rate",
          ylab = "Marginal effect of grave proximity on women as % council",
          bin.labs = FALSE) -> bin_council_covar_LIT # interaction with covariates

# > Causal mediation analysis: Political outcomes ~ genocide (via socioeconomic empowerment) ####

# >> Mediator: FHH ####
# >>> Table M-1, Model 1 ####
med.fit01 <- lm(FHEADS ~ graveDistanceCount_mean
                + mean_elev 
                + latitude 
                + dist_to_phnom_penh_km 
                + min_dist_to_provincial_capital_km 
                + dist_to_road_meters_median_km 
                + as.factor(ZONE_NAME_mode)
                ,
                data = data_commune[!is.na(data_commune$fcand_p),],
                na.action = na.omit)
out.fit01 <- lm(fcand_p ~ FHEADS + graveDistanceCount_mean
                + mean_elev 
                + latitude 
                + dist_to_phnom_penh_km 
                + min_dist_to_provincial_capital_km 
                + dist_to_road_meters_median_km 
                + as.factor(ZONE_NAME_mode)
                ,
                data = data_commune,
                na.action = na.omit)
set.seed(211)
med.out01 <- mediate(med.fit01, out.fit01, 
                     treat = "graveDistanceCount_mean",
                     control.value = 0,
                     mediator = "FHEADS",
                     robustSE = TRUE,
                     sims = 100)
summary(med.out01)

# >>> Table M-1, Model 2 ####
med.fit01b <- lm(FHEADS ~ graveDistanceCount_mean
                 + mean_elev 
                 + latitude 
                 + dist_to_phnom_penh_km 
                 + min_dist_to_provincial_capital_km 
                 + dist_to_road_meters_median_km 
                 + as.factor(ZONE_NAME_mode)
                 ,
                 data = data_commune[!is.na(data_commune$fcouncil_p),],
                 na.action = na.omit)
out.fit01b <- lm(fcouncil_p ~ FHEADS + graveDistanceCount_mean
                 + mean_elev 
                 + latitude 
                 + dist_to_phnom_penh_km 
                 + min_dist_to_provincial_capital_km 
                 + dist_to_road_meters_median_km 
                 + as.factor(ZONE_NAME_mode)
                 ,
                 data = data_commune,
                 na.action = na.omit)
set.seed(211)
med.out01b <- mediate(med.fit01b, out.fit01b, 
                      treat = "graveDistanceCount_mean",
                      control.value = 0,
                      mediator = "FHEADS",
                      robustSE = TRUE,
                      sims = 100)
summary(med.out01b)

# >> Mediator: literacy ####
# >>> Table M-1, Model 3 ####
med.fit02 <- lm(F_LIT15 ~ graveDistanceCount_mean
                + mean_elev 
                + latitude 
                + dist_to_phnom_penh_km 
                + min_dist_to_provincial_capital_km 
                + dist_to_road_meters_median_km 
                + as.factor(ZONE_NAME_mode)
                ,
                data = data_commune[!is.na(data_commune$fcand_p),],
                na.action = na.omit)
out.fit02 <- lm(fcand_p ~ F_LIT15 + graveDistanceCount_mean
                + mean_elev 
                + latitude 
                + dist_to_phnom_penh_km 
                + min_dist_to_provincial_capital_km 
                + dist_to_road_meters_median_km 
                + as.factor(ZONE_NAME_mode)
                ,
                data = data_commune,
                na.action = na.omit)
set.seed(211)
med.out02 <- mediate(med.fit02, out.fit02, 
                     treat = "graveDistanceCount_mean",
                     control.value = 0,
                     mediator = "F_LIT15",
                     robustSE = TRUE,
                     sims = 100)
summary(med.out02)

# >>> Table M-1, Model 4 ####
med.fit02b <- lm(F_LIT15 ~ graveDistanceCount_mean
                 + mean_elev 
                 + latitude 
                 + dist_to_phnom_penh_km 
                 + min_dist_to_provincial_capital_km 
                 + dist_to_road_meters_median_km 
                 + as.factor(ZONE_NAME_mode)
                 ,
                 data = data_commune[!is.na(data_commune$fcouncil_p),],
                 na.action = na.omit)
out.fit02b <- lm(fcouncil_p ~ F_LIT15 + graveDistanceCount_mean
                 + mean_elev 
                 + latitude 
                 + dist_to_phnom_penh_km 
                 + min_dist_to_provincial_capital_km 
                 + dist_to_road_meters_median_km 
                 + as.factor(ZONE_NAME_mode)
                 ,
                 data = data_commune,
                 na.action = na.omit)
set.seed(211)
med.out02b <- mediate(med.fit02b, out.fit02b, 
                      treat = "graveDistanceCount_mean",
                      control.value = 0,
                      mediator = "F_LIT15",
                      robustSE = TRUE,
                      sims = 100)
summary(med.out02b)

# >> Mediator: sex ratio ####
# >>> Table M-2, Model 1 ####
med.fit03 <- lm(SEXRATIO ~ graveDistanceCount_mean
                + mean_elev 
                + latitude 
                + dist_to_phnom_penh_km 
                + min_dist_to_provincial_capital_km 
                + dist_to_road_meters_median_km 
                + as.factor(ZONE_NAME_mode)
                ,
                data = data_commune[!is.na(data_commune$fcand_p),],
                na.action = na.omit)
out.fit03 <- lm(fcand_p ~ SEXRATIO + graveDistanceCount_mean
                + mean_elev 
                + latitude 
                + dist_to_phnom_penh_km 
                + min_dist_to_provincial_capital_km 
                + dist_to_road_meters_median_km 
                + as.factor(ZONE_NAME_mode)
                ,
                data = data_commune,
                na.action = na.omit)
set.seed(211)
med.out03 <- mediate(med.fit03, out.fit03, 
                     treat = "graveDistanceCount_mean",
                     control.value = 0,
                     mediator = "SEXRATIO",
                     robustSE = TRUE,
                     sims = 100)
summary(med.out03)

# >>> Table M-2, Model 2 ####
med.fit03b <- lm(SEXRATIO ~ graveDistanceCount_mean
                 + mean_elev 
                 + latitude 
                 + dist_to_phnom_penh_km 
                 + min_dist_to_provincial_capital_km 
                 + dist_to_road_meters_median_km 
                 + as.factor(ZONE_NAME_mode)
                 ,
                 data = data_commune[!is.na(data_commune$fcouncil_p),],
                 na.action = na.omit)
out.fit03b <- lm(fcouncil_p ~ SEXRATIO + graveDistanceCount_mean
                 + mean_elev 
                 + latitude 
                 + dist_to_phnom_penh_km 
                 + min_dist_to_provincial_capital_km 
                 + dist_to_road_meters_median_km 
                 + as.factor(ZONE_NAME_mode)
                 ,
                 data = data_commune,
                 na.action = na.omit)
set.seed(211)
med.out03b <- mediate(med.fit03b, out.fit03b, 
                      treat = "graveDistanceCount_mean",
                      control.value = 0,
                      mediator = "SEXRATIO",
                      robustSE = TRUE,
                      sims = 100)
summary(med.out03b)

# > Two-stage least squares: Political outcomes ~ x ~ genocide ####

# >> Endogenous: FHH ####
# >>> Table N-1, Model 1 ####
comm_2sls_fhh_cand <-ivreg(fcand_p ~ FHEADS + mean_elev + 
                             latitude + dist_to_phnom_penh_km + 
                             min_dist_to_provincial_capital_km + 
                             dist_to_road_meters_median_km + 
                             as.factor(ZONE_NAME_mode) | graveDistanceCount_mean + mean_elev + latitude + dist_to_phnom_penh_km + min_dist_to_provincial_capital_km + dist_to_road_meters_median_km + as.factor(ZONE_NAME_mode),
                           data = data_commune)
comm_2sls_fhh_cand_coef <- coeftest(comm_2sls_fhh_cand, vcov = vcovHC(comm_2sls_fhh_cand))

# >>> Table N-1, Model 3 ####
comm_2sls_fhh_council <-ivreg(fcouncil_p ~ FHEADS + mean_elev + 
                                latitude + dist_to_phnom_penh_km + 
                                min_dist_to_provincial_capital_km + 
                                dist_to_road_meters_median_km + 
                                as.factor(ZONE_NAME_mode) | graveDistanceCount_mean + mean_elev + latitude + dist_to_phnom_penh_km + min_dist_to_provincial_capital_km + dist_to_road_meters_median_km + as.factor(ZONE_NAME_mode),
                              data = data_commune)
comm_2sls_fhh_council_coef <- coeftest(comm_2sls_fhh_council, vcov = vcovHC(comm_2sls_fhh_council))

# >> Endogenous: literacy ####
# >>> Table N-1, Model 2 ####
comm_2sls_flit_cand <-ivreg(fcand_p ~ F_LIT15 + mean_elev + 
                              latitude + dist_to_phnom_penh_km + 
                              min_dist_to_provincial_capital_km + 
                              dist_to_road_meters_median_km + 
                              as.factor(ZONE_NAME_mode) | graveDistanceCount_mean + mean_elev + latitude + dist_to_phnom_penh_km + min_dist_to_provincial_capital_km + dist_to_road_meters_median_km + as.factor(ZONE_NAME_mode),
                            data = data_commune)
comm_2sls_flit_cand_coef <- coeftest(comm_2sls_flit_cand, vcov = vcovHC(comm_2sls_flit_cand))

# >>> Table N-1, Model 4 ####
comm_2sls_flit_council <-ivreg(fcouncil_p ~ F_LIT15 + mean_elev + 
                                 latitude + dist_to_phnom_penh_km + 
                                 min_dist_to_provincial_capital_km + 
                                 dist_to_road_meters_median_km + 
                                 as.factor(ZONE_NAME_mode) | graveDistanceCount_mean + mean_elev + latitude + dist_to_phnom_penh_km + min_dist_to_provincial_capital_km + dist_to_road_meters_median_km + as.factor(ZONE_NAME_mode),
                               data = data_commune)
comm_2sls_flit_council_coef <- coeftest(comm_2sls_flit_council, vcov = vcovHC(comm_2sls_flit_council))

# > Pre-Khmer Rouge province fixed effects ####

# >>> Table O-1 ####
data_commune %>%
  lm(fcand_p ~ graveDistanceCount_mean + as.factor(colonial_province_1920),
     data = .) -> commCand1_colonialFE
commCand1_colonialFE_coef <- coeftest(commCand1_colonialFE, vcov. = vcovHC(commCand1_colonialFE))

data_commune %>%
  lm(fcand_p ~ graveDistanceCount_mean
     + mean_elev 
     + latitude 
     + dist_to_phnom_penh 
     + min_dist_to_provincial_capital 
     + dist_to_road_meters_median 
     + as.factor(colonial_province_1920),
     data = .) -> commCand2_colonialFE
commCand2_colonialFE_coef <- coeftest(commCand2_colonialFE, vcov. = vcovHC(commCand2_colonialFE))

# >>> Table O-2 ####
data_commune %>%
  lm(fcouncil_p ~ graveDistanceCount_mean + as.factor(colonial_province_1920),
     data = .) -> commCouncil1_colonialFE
commCouncil1_colonialFE_coef <- coeftest(commCouncil1_colonialFE, vcov. = vcovHC(commCouncil1_colonialFE))

data_commune %>%
  lm(fcouncil_p ~ graveDistanceCount_mean
     + mean_elev 
     + latitude 
     + dist_to_phnom_penh 
     + min_dist_to_provincial_capital 
     + dist_to_road_meters_median 
     + as.factor(colonial_province_1920),
     data = .) -> commCouncil2_colonialFE
commCouncil2_colonialFE_coef <- coeftest(commCouncil2_colonialFE, vcov. = vcovHC(commCouncil2_colonialFE))

# >>> Table O-3, Models 1-2 ####
data_commune %>%
  lm(FHEADS ~ graveDistanceCount_mean + as.factor(colonial_province_1920),
     data = .) -> commFHEADS1_colonialFE
commFHEADS1_colonialFE_coef <- coeftest(commFHEADS1_colonialFE, vcov. = vcovHC(commFHEADS1_colonialFE))

data_commune %>%
  lm(FHEADS ~ graveDistanceCount_mean
     + mean_elev 
     + latitude 
     + dist_to_phnom_penh 
     + min_dist_to_provincial_capital 
     + dist_to_road_meters_median 
     + as.factor(colonial_province_1920),
     data = .) -> commFHEADS2_colonialFE
commFHEADS2_colonialFE_coef <- coeftest(commFHEADS2_colonialFE, vcov. = vcovHC(commFHEADS2_colonialFE))

# >>> Table O-3, Models 3-4 ####
data_commune %>%
  lm(F_LIT15 ~ graveDistanceCount_mean + as.factor(colonial_province_1920),
     data = .) -> commFLIT1_colonialFE
commFLIT1_colonialFE_coef <- coeftest(commFLIT1_colonialFE, vcov. = vcovHC(commFLIT1_colonialFE))

data_commune %>%
  lm(F_LIT15 ~ graveDistanceCount_mean
     + mean_elev 
     + latitude 
     + dist_to_phnom_penh 
     + min_dist_to_provincial_capital 
     + dist_to_road_meters_median 
     + as.factor(colonial_province_1920),
     data = .) -> commFLIT2_colonialFE
commFLIT2_colonialFE_coef <- coeftest(commFLIT2_colonialFE, vcov. = vcovHC(commFLIT2_colonialFE))

# >>> Table O-4, Models 1-2 ####

data_hh %>%
  felm(fbusinessmanager ~ graveDistanceCount
       | colonial_province_1920 | 0 | villageid, 
       data = ., 
       subset = fhh == 1,
       weights = .$hhweight[data_hh$fhh == 1]) -> biz.fhh.biv.colonialFE
data_hh %>%
  felm(fbusinessmanager ~ graveDistanceCount + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters
       | colonial_province_1920 | 0 | villageid, 
       data = ., 
       subset = fhh == 1,
       weights = .$hhweight[data_hh$fhh == 1]) -> biz.fhh.cov.colonialFE

# >>> Table O-4, Models 3-4 ####

data_hh %>%
  felm(fbusinessmanager ~ graveDistanceCount | colonial_province_1920 | 0 | villageid, 
       data = ., 
       subset = fhh == 0,
       weights = .$hhweight[data_hh$fhh == 0]) -> biz.mhh.biv.colonialFE
data_hh %>%
  felm(fbusinessmanager ~ graveDistanceCount + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters
       | colonial_province_1920 | 0 | villageid, 
       data = ., 
       subset = fhh == 0,
       weights = .$hhweight[data_hh$fhh == 0]) -> biz.mhh.cov.colonialFE

# >>> Table O-4, Models 5-6 ####

data_hh %>%
  felm(fschyrs ~ graveDistanceCount | colonial_province_1920 | 0 | villageid, 
       data = ., 
       subset = fhh == 1,
       weights = .$hhweight[data_hh$fhh == 1]) -> edu.fhh.biv.colonialFE
data_hh %>%
  felm(fschyrs ~ graveDistanceCount
       + mean_elev 
       + coords.x2 
       + dist 
       + min_dist 
       + dist_to_road_meters | colonial_province_1920 | 0 | villageid, 
       data = ., 
       subset = fhh == 1,
       weights = .$hhweight[data_hh$fhh == 1]) -> edu.fhh.cov.colonialFE

# >>> Table O-4, Models 7-8 ####

data_hh %>%
  felm(fschyrs ~ graveDistanceCount | colonial_province_1920 | 0 | villageid, 
       data = ., 
       subset = fhh == 0,
       weights = .$hhweight[data_hh$fhh == 0]) -> edu.mhh.biv.colonialFE
data_hh %>%
  felm(fschyrs ~ graveDistanceCount
       + mean_elev 
       + coords.x2 
       + dist 
       + min_dist 
       + dist_to_road_meters | colonial_province_1920 | 0 | villageid, 
       data = ., 
       subset = fhh == 0,
       weights = .$hhweight[data_hh$fhh == 0]) -> edu.mhh.cov.colonialFE

# >>> Table O-5, Models 1-3 ####

data_hh2 %>%
  felm(I(education.spending/total.spending) ~ graveDistanceCount + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | colonial_province_1920 | 0 | villageid,
       data = .,
       weights = .$hhweight) -> reg_cses9_education.spending_colonialFE
data_hh2 %>%
  felm(I(education.spending/total.spending) ~ graveDistanceCount + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | colonial_province_1920 | 0 | villageid,
       data = .,
       subset = FHH == 1,
       weights = .$hhweight[.$FHH == 1]) -> reg_cses9_education.spending_FHH_colonialFE
data_hh2 %>%
  felm(I(education.spending/total.spending) ~ graveDistanceCount + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | colonial_province_1920 | 0 | villageid,
       data = .,
       subset = FHH == 0,
       weights = .$hhweight[.$FHH == 0]) -> reg_cses9_education.spending_MHH_colonialFE

# >>> Table O-5, Models 4-6 ####

data_hh2 %>%
  felm(I(education_spending_female/total.spending) ~ graveDistanceCount + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | colonial_province_1920 | 0 | villageid,
       data = .,
       weights = .$hhweight) -> reg_cses9_education.spending_female2_colonialFE
data_hh2 %>%
  felm(I(education_spending_female/total.spending) ~ graveDistanceCount + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | colonial_province_1920 | 0 | villageid,
       data = .,
       subset = FHH == 1,
       weights = .$hhweight[.$FHH == 1]) -> reg_cses9_education.spending_female2_FHH_colonialFE
data_hh2 %>%
  felm(I(education_spending_female/total.spending) ~ graveDistanceCount + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | colonial_province_1920 | 0 | villageid,
       data = .,
       subset = FHH == 0,
       weights = .$hhweight[.$FHH == 0]) -> reg_cses9_education.spending_female2_MHH_colonialFE

# >>> Table O-5, Models 7-9 ####

data_hh2 %>%
  felm(I((education_spending_female/females_positive_spending)/1000) ~ graveDistanceCount + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | colonial_province_1920 | 0 | villageid,
       data = .,
       weights = .$hhweight) -> reg_cses9_education_spending_female_colonialFE
data_hh2 %>%
  felm(I((education_spending_female/females_positive_spending)/1000) ~ graveDistanceCount + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | colonial_province_1920 | 0 | villageid,
       data = .,
       subset = FHH == 1,
       weights = .$hhweight[.$FHH == 1]) -> reg_cses9_education.spending_female_FHH_colonialFE
data_hh2 %>%
  felm(I((education_spending_female/females_positive_spending)/1000) ~ graveDistanceCount + mean_elev + coords.x2 + dist + min_dist + dist_to_road_meters | colonial_province_1920 | 0 | villageid,
       data = .,
       subset = FHH == 0,
       weights = .$hhweight[.$FHH == 0]) -> reg_cses9_education.spending_female_MHH_colonialFE
